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6% and 30% for the computing cost in this comparison. 
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1. Introduction 

In searching for gravitational wave signals from coalescing binary compact objects, 
one commonly uses an optimal filtering technique pp. This technique consists of the 
comparison of the output signal of an interferometric gravitational waves detector with 
a family of expected theoretical waveforms, called templates. Each template depends on 
one or more parameters {Aj}. The choice of the templates in the {Aj} parameter space, 
called placement, is the purpose of this paper. We restrict ourselves to a 2D parameter 
space, considering spinless templates computed at second post-newtonian order. 

We will first describe in section |2] the motivations of our placement technique, 
comparing it with a simple uniform paving of the parameter space. Section|3]describes the 
calculation of the parameters of the parameter space portion covered by a single template. 
This portion is in our case well approximated by an ellipse. Next, section 0] treats the 
triangulation of the parameter space, a step needed by the placement, which is covered by 
section 03 Finally, performance tests are covered by section El where some real use-cases 
are considered in the context of the Virgo detector 

2. Motivations 

2.1. Portion of the parameter space covered by one template 

The comparison of a signal with one template is made through a Wiener filter [Hj: 

(a,f) =2 

This is essentially a weighted intercorrelation, a(f) being the interferometer output and 
T(f) the template. S(f) is the noise power spectral density (PSD) of the detector, / 4 and 
f s are the lower and upper limits of the detector spectral window. 

Each template is represented by a point in a multidimensional parameter space. After 
taking care of most extrinsic parameters (like time of arrival or initial orbital phase of the 
system) by maximizing the output of the optimal filter over them jj], there remain only 
two parameters, that we will call Ai and A2. Those parameters may be the masses of the 
two bodies but in general, one uses parameters derived from the masses that simplify the 
calculations. 

A template corresponding to parameters (Ai, A2) is sensitive to a signal corresponding 
to nearby parameters (Ai + SXi, A2 + 6X2)- The difference leads to a decrease in signal 
over noise ratio (SNR) with respect to the SNR obtained with a signal corresponding to 
the exact template. For an acceptable loss in SNR, each template covers a portion of 
the two dimensional parameter space. Following Owen 1 in a geometrical interpretation 
of the optimal filtering, one is able to define a distance between two templates as the 
ambiguity function maximized over extrinsic parameters, called "match" . When filtering 
a signal which has the same shape as a template of parameters (Ai + 5X\, A 2 + SX 2 ) with 
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a reference template of parameters (Ai, A2), the match is the fraction of the optimal SNR 
obtained when filtering the reference template with a signal identical in shape to itself. 

Given a minimal match MM, we can define the region of parameter space around 
a given point corresponding to a template T, the match of which, computed with any 
template corresponding to a point in the region, will be above MM. We will call the 
boundary of this region the "isomatch contour". The shape of this boundary may be 
complex, so one generally uses parameters for which it has been shown that, for high 
values of the minimal match, (MM > 0.97) the contour is closed and well approximated 
by an ellipse pQ . Throughout this paper, instead of masses, we will use chirp times r and 
defined as: 

5M _ _ -M 

T ° = 25677(7rM/ ) 8 / 3 Tl5 = 8r/(7rM/ ) 5/3 

in geometrized units (G = c = 1), where M is the total mass of the binary system, 
7] = mi7«2/M 2 is the symmetric mass ratio and fo a fiducial frequency chosen as the lower 
frequency cutoff of the detector sensitivity. Results are properly scaled to restore physical 
units. 

The calculation of the parameters of the ellipse may be done analytically for a given 
spectral density PjE]- 

The final goal of our study is to pave the parameter space with isomatch contours in 
as optimal a way as possible. This is equivalent to finding the minimal set of templates 
whose isomatch contours pave all the parameter space, without letting any hole or unpaved 
region [Zj. 

2.2. Simple paving of the parameter space 

One simple solution, already described elsewhere, is to calculate the ellipse parameters 
for the point in the parameter space where it is known to be the smallest and pave the 
space with this single ellipse [TJ, obtaining a regular tiling of the parameter space. This is 
not very different from paving a bidimensional space with circles. As was already noted 
[Zj , because of the rotational symmetry, the centers of the circles should sit at the vertices 
of regular polygons which make a regular tiling of the plane. This is only possible for 
triangles, squares or hexagons. In the first case, the centers of the circles are placed on 
the corners of an equilateral triangle, as shown in figure ^ A). It is desirable to have the 
sparsest possible circles, which means that three circles touch at one single point P. The 
surface region consisting of the points whose closest circle center is C is shown in gray. 
This is also the surface covered on average by one circle. In the triangular case, it is a 
hexagon. The set of points which belongs to this region is called the Voronoi set of C . As 
illustrated in figure [U in the case of a square tiling, the Voronoi set has a square shape 
and in the case of a hexagonal tiling, the Voronoi set has a triangular shape. It has been 
shown [Jj, as one would intuitively expect, that the most efficient tiling in the case of 
placement of circles is the triangular one. Of course, in our case, the circles are skewed 
according to the parameters of the initially calculated ellipse. 
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A) Triangular B) Square C) Hexagonal 



Figure 1. Paving of a plane with different elementary cells. The relevant Voronoi sets 
are shown in gray. 

The tiling is extended outside the parameter space to make the coverage complete. 
The ellipses, the center of which lies in a physically forbidden region (under the equal 
mass line), are shifted towards the allowed region, staying on the equal mass line, still 
ensuring the completeness of the coverage. An example is given in fig. El where the ellipse 
at the extreme right (smallest masses) represents the only computed point. 
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Figure 2. Example of a regular tiling of a parameter space. The templates are computed 
at 2 PN order in the mass range [5;50] M 0) the minimal match being 0.95, the frequency 
range [50;2000] Hz with the Virgo PSD. 



2.3. Improvements to this method 

The above simple method is very fast but, assuming that one uses the smallest possible 
ellipse, is clearly suboptimal in most cases. It gives a higher number of templates than 
would be ideally needed if one was able to calculate the shape of the isomatch contours 
at any given point of the parameter space and use those bigger shapes to cover the space. 
A second problem would then arise, since an optimal tiling of the parameter space with 
varying shapes is far from being obvious. The principle of reconstruction of exact isomatch 
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contours has been described previously jU] as well as a preliminary placement method. 

We present in this paper an extension and improvement of this method in the case 
where the elliptic approximation for isomatch contours is assumed valid. 

3. Computation of ellipse parameters 

Before doing the placement, one should be able to calculate as fast as possible the ellipse 
parameters at any given point in the parameter space. This is done by 

• Calculating the ellipses at a chosen set of points (we obtain "seed ellipses"). 

• Triangulating the parameter space with this set. Actually, as we will see, those two 
steps are closely linked. We give in appendix a short tutorial about triangulation and 
computational geometry. 

• Interpolate linearly ellipses at any point using the previously calculated seed ellipses. 
This step is much faster than an analytical computation. 

3.1. Computation of seed ellipses 

The seed ellipses are computed using the algorithm included inside the LIGO Analysis 
Library (LAL) [13J. This algorithm uses the procedure described in jH]. The metric 
components used to find the parameters of the ellipse are calculated using the moments 
of the PSD curve. 

3.2. Triangulation and interpolation 

The triangulation of the parameter space deserves hereafter a section by itself. Once it 
is computed, each point P in the parameter space belongs to one and only one triangle 
whose corners are three seed points. One is able to interpolate linearly the shapes (resp. 
metric parameters) of the three seed ellipses to obtain the parameters of the ellipse (resp. 
metric) at point P (see fig. Ej). 

4. Triangulation of the parameter space 

The triangulation of the parameter space is done using standard techniques known in 
computational geometry. The notions necessary to understand the present study are 
explained in appendix. The base algorithm used is known as the Bowyer- Watson |14j|15j 
algorithm. 

4-1. Triangulation algorithm adapted to the CB parameter space 

The Bowyer- Watson algorithm is quite simple but needs adaptation to our problem. We 
need to take care of the fact that the borders of the parameter space are not convex and 
we need to choose which points to use for the triangulation. 



Variable placement of templates technique for binary inspiral searches 



6 




Figure 3. Linear interpolation between analytically calculated ellipses (E%, E2, E3), 
that may differ in size and orientation, at an arbitrary point P in a triangle. 



The main idea of our adapted algorithm is to start from an existing triangle at the 
corners of which sit three already calculated ellipses {Ei, E 2 , E 3 } and subdivide it only 
if necessary, i.e. if for any point P inside the triangle, the ellipse linearly interpolated 
between {E\, E2, E3} is different enough from the one calculated using the metric at that 
point. Let Ei be the interpolated ellipse and E c the calculated one. a being the measure 
of the surface of a out the surface of Ei that does not intersect E c (fig. HJ), the variable 
describing the difference between Ei and E c has been chosen as the proportion 



It was not deemed necessary to also take into account the surface of E c that does not 
intersect E iy because if a out is null, the interpolated ellipse is completely inscribed inside 
the calculated one and we are simply going to make a more dense placement at a later 
stage. A limit is set on this variable to stop the subdivision of triangles. 




Figure 4. The variable indicating the difference between the two interpolated and 
calculated ellipses at the same point is the proportion of the surface of the interpolated 
ellipse not common to both ellipses. 
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4-1.1. Division of an existing triangle Given an existing triangle, a choice has to be 
made on the points appropriate for its subdivision. Ideally, one would use the points 
which have the highest proportion p. It is however impractical, and very expensive in 
terms of computing power to test all the points in a triangle to find the one with the 
higher p. We chose to test only the middle points of each segment forming the triangle. 

Each of these three points is inserted and used to subdivide the triangle following 
a Delaunay method, but considering only the triangle, not the adjacent ones that may 
exist in the ongoing triangulation process. If the middle point of a segment is outside 
the parameter space, it is replaced by the closest point on the border, perpendicularly to 
the segment (fig. EJ). Some peculiar situations (two middle segment points outside the 
parameter space for example) are taken into account. All subtriangles generated outside 
the parameter space are removed. 




Figure 5. Replacement of a point M outside the parameter space by a point M' on 
the border close by (a). A peculiar case, which is taken into account, is also shown (b) 



From the description above, it is obvious that the final triangulation will not be 
strictly speaking a Delaunay one, since we use the Delaunay criterion only locally for a 
triangle subdivision. 

4-1.2. Global algorithm view We start from the triangle formed in the (tq, T1.5) parameter 
space by the three angular points corresponding to (m min ,m min ), (m min ,m max ) and 
(rn m ax,rn max ), where m m i n and m max are respectively the minimal and maximal masses 
of the binary system members considered. 

The triangle is then recursively subdivided as described above. Since there is a limit 
on the p proportion of each inserted point, the subdivision will stop naturally when the 
mesh becomes dense enough. These successive refinement steps are illustrated in figure 
IH1 In order to avoid too big a number of calculated ellipses, and to limit the computing 
time, the number of refinement steps has been arbitrarily limited to 7. It was verified that 
this doesn't bring any problems, except in the lower left corner of the parameter space, 
corresponding to high masses (above 10 M ) for both objects. In that case, the placement 
may be somewhat wrong but a posteriori Monte-Carlo tests show an undercoverage not 
exceeding 2% of the total parameter space surface. 



Variable placement of templates technique for binary inspiral searches 



8 




Refinement step 1 Refinement step 5 



Figure 6. Mesh refinement steps. Starting from a triangle enclosing the parameter 
space, insert points and retriangulate while the proportion of discrepancy p between 
interpolated and calculated ellipse is greater than a limit pu m - 

As may be noted on figure EH the tesselation of the parameter space is extended 
in the physically allowed region to avoid some extrapolation side-effects in the following 
placement procedure. 

4-2. Extrapolation outside border of the parameter space 

Each ellipse calculated for the placement procedure described hereafter is actually 
interpolated inside one of the triangles found during the triangulation step. If the point 
considered by the placement is outside of the tesselated (triangulated) part, it doesn't 
belong to any triangle a priori. We will see that the placement procedure needs to spill 
over the strict borders of the parameter space to ensure complete coverage, and it may 
happen that a determination of ellipse parameters is needed outside the tesselated part. 
Furthermore, the calculation of the metric is impossible in the disallowed (physically 
forbidden) region under the equal mass line in the (r , T1.5) parameter space. Therefore, 
we cannot triangulate that region since we cannot calculate true ellipses or contours in it. 

Thus, we need to provide a way to extrapolate the ellipse parameters outside the 
strictly tesselated part of the space. As will be seen later, the final step of the placement 
procedure consists of shifting the points found in the forbidden region so that they fall 
in the allowed one. But extrapolation is needed all around the space border during the 
placement, albeit in a limited area. 

For a given point A outside the parameter space, it is natural to associate it with 
the closest triangle of the tesselation. The word "closest" should be taken with care, as 
closest in euclidian distance doesn't mean more adequate for our purposes. We consider 
only the triangles which border the tesselated region, i.e. which have one side that is not 
common to another triangle, thus delimits the border of the tesselation (fig. EJ). Once the 
triangle associated with the point A is determined, one can do a linear extrapolation of 
the ellipse parameters, as for the points inside the triangle. 
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Figure 7. Association between a point outside the triangulated part of the parameter 
space and a triangle on the border 

The choice of the triangle associated with a given point A is done as follows. We 
define the vectors (3ij which join two successive vertices Bi and Bj lying on the border 
of the tesselated part of the parameter space. Each vertex Bi is associated with a vector 
ki whose direction is pointing towards the outside of the space and is an average of the 
normal to two consecutive vectors fiki and fy. 

A point A will be associated with the triangle containing the vertices Bi and Bj if it 
is located in the domain delimited by and the two lines defined by (Bi, ki) and (Bj, kj). 
An example of such a domain is shown in gray in figure 

Clearly, the very simple extrapolation we describe is valid only for the points close to 
the space boundary. The lines (Bi, ki) will cross and it is not possible to associate a point 
and a triangle beyond those crossings. Furthermore, the validity of the extrapolation is 
not guaranteed for points pushed away from the boundary of the parameter space. In our 
case, where we marginally extend the calculation of the metric outside the borders, this 
shows not to be a problem. 

4-3. Results in concrete cases 

Figure |H1 shows triangulation in a few real cases. 

• n^min and m max are the minimal and maximal masses of the parameter space 

• Fi and the lower and higher frequency cutoffs used for the generation of templates 

• PN is the order of the post-newtonian expansion 

• NStep max is the limit imposed on the number of triangulation steps 

• N Step final is the number of steps effectively needed to satisfy the surface proportion 
condition for all the ellipses generated, without reaching the NStep max limit 

• Nt is the number of calculated points in the triangulation to reach the NStep max or 
N Step f i na i limit 

• The noise spectral density used was a Virgo-like one, shown on fig. 
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Figure 8. Examples of triangulations representing real use cases for CB searches in 
Virgo. The black line represents the border of the parameter space. The triangulation 
area is extended outside to prevent extrapolation problems. See the text for an 
explanation of computing conditions. 




1 10 10* 10 3 10 4 

Hz 



Figure 9. Noise spectral density for the real use cases. 



5. Placement 

5.1. I somatch properties 

Once the triangulation and seed ellipses have been generated, the placement is done in 
two stages. Both rely on basic properties of isomatch contours described in j^j, namely: 

• The match symmetry between two contours. If T\ and T 2 are two normalized 
templates, one has (Ti,T 2 ) = (T 2 ,Ti) = M. Thus, the point in the parameter space 
corresponding to T\ is located on the isomatch contour of value M corresponding to 
T 2 , and conversely, the point corresponding to T 2 is located on the isomatch contour 
of value M corresponding to T\. In practical computations, the match symmetry may 
not be absolutely verified because in general one maximizes over the initial phase of 
one template (say Ti), which is not done for the signal (T 2 ). This has proven to 
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be negligible for smooth variations of the metric throughout the parameter space, 
which is roughly the case in our tests using the LAL, except perhaps for high masses, 
> 10 M . 

• To place an ellipse with respect to another in an optimal way, one introduces a guiding 
ellipse. This allows to place three ellipse sets (fig. [TUJ). The three ellipses intersect 
at the center of the guiding ellipse. 

In the course of the running of the algorithm, if two of the ellipses are placed, the 
third one may be positioned naturally on the border of the guiding ellipse by maximizing 
the surface of the triangle formed by the centers of the three ellipses. 




Figure 10. Placement of 
three contours in the {tq,ti^) 
plane 



Figure 11. First stage 
of placement along the equal 
mass line 



5.2. First stage 

The first stage consists of a side by side placement along the equal mass line starting from 
the (m max , m max ) point (fig. ITTj) . Unlike in the simple placement case where this was 
avoided, it is the most efficient way of paving while only one ellipse is needed to cover 
the parameter space along the direction of the semi-major axis of the ellipse, an almost 
vertical direction in most of our cases. 

The principle is described in figure El Starting from an ellipse Ei the center of which 
Ci lies on the equal mass line, a choice is made (explained hereafter) of the position C g 
of the center of a guiding ellipse along the border of E^. Because of the isomatch contour 
properties stated above, Cj lies on the guiding ellipse. It is also on the equal mass line. 
Ci + \ is the other intersection of the guiding ellipse and the equal mass line. 

The position C g is chosen between Ai and A2 limits on the Ei ellipse, in such a 
way that the surface of the {Cj, CV^C^} triangle is maximized. C3 is the location 
of the center of a potential ellipse E3 that would form with Ei and the next ellipse 
E i+ i a three ellipse set optimally placed (with the placement conditions imposed by the 
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Choice of the position of the guiding ellipse Position of the next ellipse 

Figure 12. principle of the first stage of placement along the equal mass line 

parameter space lower boundary correponding to the equal mass line). The Ai and A2 
limits are chosen empirically and are subject to the influence of numerical errors as well 
as interpolation/extrapolation errors. 

The next ellipse E i+ i is then placed at position C i+ i. E i+ i and E { should ideally 
intersect at two points Sh and Si, Sh being equal to the center of the guiding ellipse C g and 
Si being in the physically forbidden region underneath the equal mass line. Because of the 
curvature and variation of the metric, it may happen that E i+ i and Ei do not intersect. 
In that case, the position C i+ i of E i+1 is shifted towards Cj along the equal mass line 
until the point C g comes on E i+ i. 

The first stage placement algorithm stops when Sh falls inside the parameter space, 
which means that two ellipses are needed to cover the parameter space in the vertical 
direction. 

5.3. Second stage 

The second stage of placement consists of the coverage of the parameter space line by 
line, as was described in 0. 

• One starts from a three ellipse set placed optimally at a point D . 

• Then place iteratively ellipses using successive guiding ellipses that follow rules 
defined in section 15.11 The placement is done alternatively on the left and on the 
right of the line of guiding ellipses, and successively above and below the initial point 
D . 

• One obtains a two-line set crossing the parameter space (fig. ITHJ) . Among the external 
crossings of the generated ellipses of one of the lines (called 71$), a point D\ is chosen 
and the process is iterated. 
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• At each step, only one of the lines is kept, the other being approximately 
superimposed with a line generated at the previous step (fig. IT5]) . 

• The starting point of each two-line set Di for the step j is chosen among the 

as the point outside the parameter space and not in the physically forbidden region 
which is the closest to the border of the parameter space. Other choices have led to 
the observation of variations in the direction of two successive lines, giving holes in 
the coverage of the space. 




Figure 13. Independent placement of the third line starting from a crossing point on 
the edge of a first two line set. 

The starting point for the first line building of the second stage is the first intersection 
point 5h found in the first stage that is inside the parameter space. The placement ends 
when no ellipse from a line covers any part of the parameter space. 

5.4- Correcting points felt outside of the parameter space 

Once the first two stages are finished, a cleaning is performed to remove superfluous 
ellipses that do not cover any part of the parameter space. 

It is not possible to do it beforehand because it is not obvious if a given ellipse 
covers or not a part of the parameter space before it is actually placed. Its center may lie 
outside of the parameter space but a small part of the ellipse may still cover a portion of 
the parameter space. 

A position correction is also done on ellipses which, while covering a portion of the 
parameter space, have their center in the physically forbidden region. Those ellipses are 
shifted following the guiding contour used for their generation until they fall on the equal 
mass border. 

5. 5. Examples of computed placements 

Figures ITU and IT5l show a few real use-cases of placement. 

• mmin and m max are the minimal and maximal masses of the parameter space 
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m min = 1M , m m ax = 30M©, MM = 0.95, F t = 40 Hz, F h = 2000 Hz.PiV = 2, N T = 152, N P = 11369 



Figure 14. Example of placement obtained for real world CB searches in Virgo. The 
black line represents the border of the parameter space. Three portions of space are 
shown and color of ellipses is varied to help viewing the shapes. See the text for an 
explanation of computing conditions. 




Fi = 40 Hz, F h = 2000 Hz, N P = 1011 F t = 40 Hz, F h = 100 Hz, N P = 226 

Only first stage needed to pave all space 
m mm = 1.95M , m max = 4.O5M , MM = 0.95, PN = 2 



Figure 15. Example of placement obtained for real world CB searches in Virgo. The 
black line represents the border of the parameter space. Color of ellipses is varied to help 
viewing the shapes. See the text for an explanation of computing conditions. 

• MM the minimal match, Fi and Fh the lower and higher frequency cutoffs used for 
the generation of templates 

• PN is the order of the post-newtonian expansion 

• Nt is the number of calculated points in the triangulation 

• Np is the number of points found in the placement. 

• The PSD used was a Virgo-like one (fig. EJ) 
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6. Performance tests 



6.1. Number of templates with the simple placement algorithm 

The number of templates needed for complete space coverage represents a simple 
performance estimator. An estimation of this number was already given jH] by computing 
the ratio between the volume of the parameter space and the proper volume covered by 
a single template. It was supposed that the packing algorithm used was a square (or 
hypercubic in D dimensions) one. The proper volume is then, in 2 dimensions 

AV = 2(1 - MM) (4) 

and in the triangular lattice case (hexagonal Voronoi sets), which was used in our simple 
algorithm for D — 2 

3\/3 



AV 



[1 - MM) 



(5) 



Table^shows the numbers we found by using an actual Virgo noise power spectral density. 
V ps being the volume of the parameter space, N p quare is the reference number computed 
as in jH] assuming a square packing, ]Sf r p f Trian9 assumes a triangular packing algorithm 
and ]\f- p imple is the actual number found with our simple algorithm described in paragraph 
2.21 which also produces a triangular lattice. Edge effects appear clearly, as the smallest 
the volume of the parameter space, the largest the difference between ]\f p e f Trmn 9 anc i 



N- 



simple 



m min (M R ) 


v ps 


j^ref Square 


jyrefTrtang 


simple 


0.5 


1.43 x 10 4 


1.43 x 10 5 


1.10 x 10 5 


1.28 x 10 5 


1 


2.57 x 10 3 


2.57 x 10 4 


1.98 x 10 4 


2.59 x 10 4 


3 


1.38 x 10 2 


1.38 x 10 3 


1.06 x 10 3 


1.80 x 10 3 



Table 1. Comparison of the number of templates obtained with a simple algorithm 
(triangular packing) and a theoretical number (ratio between parameter space volume 
and proper volume asuming a square or triangular packing). The minimum mass 
varies from 0.5 M to 3 M Q . Other conditions : m max = 30M Q , = 30Hz, 

f rnax = 2000Hz, MM = 0.95, Virgo-like PSD. 



6.2. Performance gain 

With the grid of templates coming out of the new placement algorithm, one can expect a 
gain in the total computational cost needed to perform a search over the defined parameter 
space with respect to the simple placement algorithm (paragraph 12. 2|) . This gain is not 
easily quantified because it depends on the specific search algorithm and on aspects that 
do not depend on the computational algorithm itself, such as I/O. But it may be estimated 
in at least two ways: 
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• firstly, by the gain in the overall number of templates coming out of the placement 
algorithms (method A). 

• secondly, by modeling the "standard" method for doing the optimal filtering and 
searching for an approximation of the gain (method B). 

The optimal filtering technique and an estimation of the computational cost are 
described by Schutz in [HI]. An approximation of the cost (number of floating point 
operations) for analyzing a set of N to t data values for a given template of length N s , and 
a fractional overlap x of successive data set chunks is: 

Nfi ap = — + (6) 
J F 1 — X X 

A discussion on the optimal value of x is made in [16] . but it does not take into 
account the I/O costs, as well as exchanges of data between memory and processor, 
which is found to be critical in our case. Therefore, as explained in |17j . we choose x so 
as to roughly optimize the length and the number of the vectors to be exchanged between 
the core memory and the CPU. Starting from the expression above and fixing x = \ for 
each template, it may be shown that the approximate total number of operations needed 
to analyze a set of M templates is given by: 

M 4 

Nj lop = 6N tot J2[ln2(2f s T i ) + -] (7) 
i=i ' } 

where f s is the sampling frequency and 73 the length of an individual template. This 
leads to consider a computing performance estimator of the form 

4, 



f = ^[/n 2 (2/ s r,) + 



8=1 

Since we only want an approximate expression, we consider n = To(i), where To(i) is 
the newtonian chirp time of the coalescing binary producing a given template. We made 
comparisons between the placement produced by the simple method of paragraph 12 . 21 and 
the full placement method. Tables [21 El an d El show the gain on the number of templates 
and the gain on the performance estimator £. The conditions of the tests were varied but 
the base conditions were the following: 

• minimal mass m min = 1M , 

• maximal mass m max = 3OM 

• minimal frequency for template generation / TO j„ = 30Hz 

• maximal frequency for template generation f m ax = 2000Hz 

• minimal match MM = 0.95 

• power spectral density close to the final Virgo one (fig. EJ) 
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In the tables, N r p quare is the reference number of templates, computed as in 
jH] assuming a square packing, as explained above in section 16.11 N P represents the 
number of templates found by the placement to cover the parameter space, Gn p = 
^.simple _ full simple - g ^ e g a j n j n the number of templates obtained when going from 

the simple placement to the full placement method, Nt is the number of seed templates 
necessary to triangulate the parameter space, is the gain in performance estimator. 
Unless otherwise noted, the triangulation process was stopped after 7 steps of refinement, 
which was shown in the section 14.1.21 not to bring problems. 



MM 


jy-ref Square 


^simple. 


N f p uU 


GjSTp 


Nt 


^simple 






0.90 


12840 


14047 


10746 


23.5% 


395 


2.46 x 10 5 


1.91 X 10 5 


22.4% 


0.95 


25680 


26183 


20161 


23.0% 


381 


4.58 x 10 5 


3.57 X 10 5 


21.9% 


0.98 


64200 


61144 


47183 


22.8% 


394 


1.06 x 10 6 


8.34 X 10 5 


21.6% 



Table 2. Gain on the number of templates and on performance coefficient with respect 
to a simple algorithm and varying the minimal match. Other conditions: m m j„ = 1M , 
m max = 30M Q , f mm = 30Hz, f max = 2000Hz, Virgo-like PSD. 



In table |2J the minimal match was varied from MM = 0.90 to MM = 0.98, keeping 
the other parameters equal. As can be seen, an average performance gain of roughly 
22% is achieved. It may be noted that the number of templates may also be used as a 
performance estimator, giving numbers very similar to £. 



m min (M ) 


jy simple 


Np UU 




N T 


^simple 




G, 


0.5 


129507 


113531 


12.3% 


335 


1.28 x 10 6 


1.13 x 10 6 


11.9% 


1 


26183 


20161 


23.0% 


381 


4.58 x 10 5 


3.57 x 10 5 


21.9% 


3 


1829 


1106 


39.5% 


416 


1.65 x 10 4 


1.01 x 10 4 


38.7% 



Table 3. Gain on the number of templates and on performance coefficient with respect to 
a simple algorithm. The minimum mass varies from 0.5 M to 3 M . Other conditions: 
m rnax = 3OM , f min = 30Hz, fmax = 2000Hz, MM = 0.95, Virgo-like PSD. 



In table El only the minimal mass of the stars, hence the size of the parameter space, 
was varied from m min = 0.5M Q to m min = 3M Q . The gain is naively expected to increase 
with the size of the parameter space. The bigger the parameter space, the higher the 
variation of metric, hence the bigger the variation in size of the ellipses. The results 
shown in table 01 vary in the opposite direction. This is explained by edge effects, where 
the influence of ellipses covering a small part of the parameter space, on or outside the 
border, and the way they are placed, play a dominant role. 

Finally, table |U shows the results for a variation in the frequency range. The mass 
range was limited to [1; 5] M because for high masses we are reaching the limits of the 
numerical relevance of the metric calculation. 
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fmin fmax (Hz) 


-^simple 


N fuu 


Gn p 


N T 


^simple 






30-100 


2191 


1948 


11.1% 


90 


4.09 x 10 4 


1.37 x 10 4 


11.2% 


100-2000 


986 


833 


15.5% 


167 


1.37 x 10 3 


1.17 x 10 3 


15.2% 


30-2000 


12641 


11820 


6.6% 


61 


2.33 x 10 5 


2.18 x 10 5 


6.5% 



Table 4. Gain on the number of templates and on performance coefficient with respect to 
a simple algorithm. Two limited frequency ranges are tested: [30; 100] Hz and [100; 2000] 
Hz. The mass range is [1;5] M Q . Other conditions: MM = 0.95, Virgo-like PSD. 



It may be noted that in practical algorithms, templates will be grouped by groups of 
similar length. The expression of £ (equation |HJ) will take a linear form as a function of 
the number of templates. This should bring the gains we obtained for the performance 
estimator closer to the ones obtained with the number of templates. 

6.3. Coverage tests 

The metric calculation is approximate, especially in the high mass region, where there is 
yet no good model of coalescence. It is therefore important to do independent tests on 
the covering efficiency Monte-Carlo tests were performed by testing randomly scattered 
points over the parameter space. The distribution of position is uniform in (r , 71.5) 
parameters. For each point, the corresponding waveform is computed and the match with 
the templates of the bank is calculated, retaining the highest. Actually, only the subset 
of templates which are closer than a given distance to the point, in the metric sense, are 
considered. 

The chosen conditions in terms of masses, frequencies and minimal match are the 
standard ones described in 16.21 Figure ITT)! A shows the distribution of test points over the 
parameter space, while figure ITTJlB shows the distribution of points the match of which is 
lower than the specified match (0.95 in our case). 

The low match points (with match M < 0.95) represent 1.6% of all the test points. 
There are two possible reasons for the presence of these points. The first is the presence 
of holes in between ellipses, due to suboptimal placement, the second is a possible 
miscalculation of the metric in some peculiar cases, for example for high mass binaries. 
Finer Monte-Carlo tests were performed in small regions relevant for the two cases, and 
low match point positions were superimposed with isomatch ellipses. The first case is 
illustrated with figure El where it is clearly seen that most of the low match points fall 
in existing holes of the placement. 

The second case is illustrated in figure fTSl The test is made with points chosen in the 
region r G [6.8; 7.1] and n.5 G [1.7; 2] (region named cthm), the points being inside the 
parameter space. It is clear from the picture that all the points of o~hm fall well inside an 
existing ellipse, hence they should have had match M > 0.95 if the metric was correctly 
calculated. This situation is explained by the miscalculation of the ellipse orientation, as 
is illustrated in figure IT9l In this figure, the points with M > 0.99 were considered, and 
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all test points | Entries 10O0FI test points M<0.95 | Entries 16T| 




Figure 16. Monte-Carlo distribution of points in coverage tests. A: all the points, B: 
points with match M < 0.95 with the closest template. The conditions are the standard 
ones. 




Figure 17. Superposition of a small part of calculated isomatch ellipses in placement 
and test points with match M < 0.95. The test points clearly fall in holes formed by 
locally incorrect placement. 

they form a figure clearly showing the wrong orientation of the computed ellipses (several 
colors depending on the value of the match were used, the darker the points the higher 
the match). 

Figure EH compares the distribution of the test points match for the full placement 
algorithm and for the simple placement algorithm. The simple placement is clearly 
suboptimal, but ensures a complete covering of the parameter space while the optimality is 
better for the full algorithm, though it does not cover perfectly the parameter space, at the 
level of a few percent undercoverage. Figure illustrates the influence of miscalculation 
of the metric. Superimposed to the distribution of the match in the full placement case, is 
the distribution of the match for (Thm- This distribution was scaled down proportionately 
to the surface of o~hm region versus the surface of the parameter space to show its 
contribution to the overall distribution. 

In general, the two effects, miscalculation and misplacement, are both present with 
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Figure 18. Superposition 
of a small part of calculated 
isomatch ellipses in placement 
and test points with match 
M < 0.95. The gray box 
surrounds the test region. 



Figure 19. Superposition 
of a small part of calculated 
isomatch ellipses in placement 
and test points with match 
M > 0.99. This illustrates the 
wrongly calculated orientation 
of the ellipses, leading to a 
wrong placement. 
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Figure 20. Comparison of 
the match distribution for the 
simple and full placement. 
The red vertical line shows the 
requested minimal match, 0.95 
in this case. 



Figure 21. Estimation of the 
contribution of points in a re- 
gion where ellipses correspond 
to high mass with miscalcu- 
lated metric. The red vertical 
line shows the requested mini- 
mal match, 0.95 in this case. 



various strength throughout the whole parameter space. Miscalculation is due to wrong 
approximation of the metric and/or approximations in the triangulation and interpolation 
steps of the placement algorithm. 

In figure ^JB, a clear accumulation of low match points seems to occur in the high 
mass region. To confirm this, a test was made with a mass range [m min ; m max ] = [1; 1O]M . 
The match distribution for this test, superimposed on the match distribution of the entire 
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parameter space ([m min ; m max ] = [1;3O]M ) and on the match distribution of cthm, is 
shown in figure 1221 It is clear from this figure that the main source of low match points 
is the miscalculation of the metric for high masses. 




match 

Figure 22. Comparison of match distributions for [m m i n ;m max ] — [1;30]M Q and 
[m m i„;m raal ] = [1;10]M Q . Also included is the distribution corresponding to a 
specifically high mass region in the parameter space, scaled proportionately to its surface 
with respect to that of the whole parameter space. 

In order to get an idea of how to easily overcome these problems, one can calculate the 
proportion of bad match test points as a function of a varying minimal match value M , for 
a given placement. This corresponds to enlarging the ellipses obtained with a placement 
with an initial minimal match 5W - Figure 1221 shows the variation of the proportion p of 
test points with match lower than M versus M . From this figure, given a desired bad 
match points proportion, one gets an estimation of the effective minimal match reached. 

A question may be raised about the robustness of the algorithm, i.e. is the algorithm 
adequate for real, noisy data. It is very difficult to assess an "absolute" robustness of the 
algorithm, because of three main points : 

• the difference between the true contour and the calculated ellipse, especially for low 
matches, which may in some circumstances push the algorithm to its limits. 

• the fact that the algorithm is not robust in the case of large and fast variations of 
the metric. 

• the difference between calculated and interpolated ellipses that may, though not in 
large proportions, affect the algorithm. 

There is a need for tools that run online and verify the relevance of the computed grid 
banks. In the case of problems, it is always possible to switch to the simple algorithm. 
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Figure 23. Proportion of bad match test points as a function of an effective minimal 
match value for a given placement. Only statistical errors are reported. 



6.4- Speed tests and recomputation of the placement 

All the tests were performed on a Linux 2.4 GHz Pentium IV workstation and we present 
in table 03 the computation time in seconds needed for each placement, in increasing 
number of generated grid points. The time is divided in two, corresponding to the two 
main steps of the algorithm, namely first the triangulation and generation of seed contours 
and second the placement itself. There is a rough proportionality between the number of 
final grid points and the time, with a quasi constant term corresponding to the first step 
(the number of generated seed contours being always of the same order of magnitude). 





MM 


N T 


Np UU 


Time (s) 


3-30 


0.90 


418 


643 


60 


3-30 


0.95 


416 


1106 


67 


3-30 


0.98 


402 


2362 


84 


1-30 


0.90 


395 


10746 


203 


1-30 


0.95 


381 


20161 


294 


1-30 


0.98 


394 


47183 


575 


0.5-30 


0.90 


335 


59103 


742 


0.5-30 


0.95 


346 


113531 


1287 



Table 5. Computation time for different placements on a 2.4 GHz Pentium IV Linux 
"workstation. 



The frequency of recomputation of the placement is still under consideration in Virgo. 
It depends on the change rate of the shape of the sensitivity curve over time, the stability 
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of which is not yet fully assessed for future science runs. The numbers given in table |5] 
may seem too large for a frequent recomputation, for instance every 15 minutes, in the 
case of large volume parameter spaces. Though such a frequency is not expected for the 
final Virgo science runs, we may need to consider a parallelization of the algorithm. The 
part of the algorithm that could be parallelized efficiently is the placement part, but one 
should not expect more than an estimated factor 2 to 5 improvement in overall computing 
time, due to the sequential nature of the algorithm. Indeed, in one line of ellipses, ellipse 
number n may not be placed before ellipse number n — 1. Only the placement of complete 
lines may be somewhat decorrelated. 

7. Comparison with previous studies and perspectives 

Beside very important pioneering efforts ppjE] the results of which are now widely used, 
several previous studies were done for the template placement problem. We believe that 
our method is somewhat complementary to them. For example, the placement algorithm 
used in [IE] for extended hierarchical searches is based on a square tiling. This is justified 
in this case by the low minimal match value used (r = 0.8), which gives very irregularly 
shaped contours. Our method could probably be adapted to such a case by applying 
methods such as in 5: to determine the shape of the contours, but an important effort 
has to be made to improve the speed of the shape reconstruction algorithm, which is going 
to be one of the main limiting factors. 

Another example is the paper of Arnaud et al. where authors devise a 2D tiling 
method and test it in the case of supernova ringdown signals. It is very difficult to make 
a direct comparison between this algorithm and ours. The very large parameter space 
curvature described by Arnaud et al is likely to bring some holes if we apply directly our 
tiling method to ringdown signals. This would imply the need for an improvement to 
our placement procedure. On the other hand, the Arnaud 2D tiling method was not yet 
applied to the case of a (70,71.5) inspiral parameter space and it is not clear what would 
be the result in terms of speed and possible overcoverage. 

The computational geometry tools that we used are still valid in higher dimensional 
spaces. It may be tempting to consider the extension of our algorithm to multidimensional 
searches. In that case, the main challenge would be to improve the algorithm speed, since 
the number of contours in nD is roughly going as 

N n « N 2 2 (9) 

Where N 2 is the number of contours obtained in 2D. This is of course a "worst case" 
scenario where the granularity is the same (and high) in all dimensions. 

8. Conclusion 

We presented a technique for doing the placement of isomatch ellipses on a template 
parameter space using triangulation and interpolation of seed ellipses. A comparison is 
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done with a simple regular triangular tiling using a single ellipse. This comparison shows 
an improvement between 6% and 30% depending on the mass range and frequency range. 
Some coverage tests were also performed that show a few percent undercoverage of the 
parameter space, mainly in the high mass region. This undercoverage seems to come from 
the miscalculation of the metric for high masses. Finally, speed tests were made. 
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Appendix: A few notions of computational geometry 

Since computational geometry is not very commonly used in our field, we will give a very 
short introduction to the notions useful for the present study It is in no way exhaustive 
or pretending to be accurate. More details may be found in [XQ| or [TTj . 

Definition of a triangulation 

Given a set S of points in a euclidian space, 2-dimensional in our case, we would like 
to subdivide the space into a set of triangles, each triangle being formed by three points 
from S. Any point P in the space belongs to (is included into) one and only one triangle. 
This is however not enough and the properties of the set of triangles should be the ones 
of a triangulation. 

Some definitions first. Let's consider a set of (n + 1) affinely independent points in 
an n-dimensional euclidian space M. n . 

• The convex hull of a set of points is the minimal convex set containing all the points 
(imagine a rubber band stretched so that it encompasses all the points). 

• A simplex is the convex hull of a set of n + 1 points (a line segment in ID, a triangle 
in 2D, a tetrahedron in 3D,...). 

A triangulation T of the set of points S in M. n is a subdivision of M" into n-dimensional 
simplices such that: 

• The set of points that are vertices of the simplices coincides with S. 

• Any two simplices in T intersect in a common face, only one vertex or not at all. 

• The convex hull of S defines a domain Q in M n . If K is a simplex, then 



We illustrate the above definition by showing what is and what is not a triangulation in 
a 2-dimensional space in figure |2U for a given set of points. 

Voronoi diagram 

A triangulation is not unique, as may be seen in figure "2"3 All triangulations are not 
equivalent for a given problem. There is a need to define a criterion of suitability. The 
most commonly used criterion is the Delaunay criterion which constraints the compactness 
of the triangles and will be explained later. It is linked to the so called Voronoi diagram. 
Given S a set of points Pi in a <i-dimensional space, the Voronoi diagram is the set of cells 
Vj associated with each point Pj and defined as 



Where d is the euclidian distance between two points. In other words, V« is the locus of 
points in W 1 closer to Pi than to any other point of S. It has been shown 1*2] that the 
geometrical dual of the Voronoi diagram is a triangulation, the Delaunay triangulation 




(10) 



Vi = {Pe R n such that d(P, Pi) < d(P, Pj), Vj ^ i} 




(fig. EED- 
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Figure 25. Two examples of valid triangulations for the same set of points. Intuitively, 
the a) case is "better" than the b) case 




Figure 26. a) the Vorono'i diagram of a set of points and b) the dual of the Voronoi 
diagram, the Delaunay triangulation 



Delaunay triangulation 

The Delaunay criterion states that the open circumdisk (in 2 dimensions, circumsphere 
in n dimensions) of a triangle (simplex) contains no point from the set. The example in 
figure EH shows a triangulation not satisfying the Delaunay criterion. Among all possible 
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a) b) 

Figure 27. Example of a) a triangulation not satisfying the Delaunay criterion (the 
point P is inside the circumdisk of K) b) satisfying it. A few circumcircles of triangles 
are shown 

triangulations, the Delaunay triangulation 

• maximizes the minimum angle formed by the faces of the triangles 

• minimizes the biggest diameter of the circumcircles associated with the triangles 

Intuitively, this would mean that the Delaunay triangulation produces the more 
"compact" triangles. 

A simple algorithm 

Based on the previous definition of the Delaunay criterion, it is possible to devise a simple 
algorithm to compute a triangulation based on a set of points. It is called an incremental 
algorithm, or Bowyer- Watson algorithm |14j|15j . 

The algorithm is incremental in the sense that the points of the set S are added one 
by one, recomputing a triangulation at each step. The process starts by the generation of 
a supertriangle that encompasses all the points in S. At the end, all triangles that share 
one edge with the supertriangle are removed. The addition of one point is illustrated in 
figure |2i 

To add one point P, all the triangles whose circumcircle contains P are first removed. 
The resulting hole in the triangulation has a polygonal shape. New triangles are formed 
between P and the outside edges of the polygon. 
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Figure 28. Incremental algorithm, a) add a point P to an existing triangulation, b) 
remove all triangles whose circumcircle contains P, c) obtain a polygon enclosing P, d) 
triangulate only this polygon 



